function [select_windows_th_act,select_windows_deriv_act,windows_TP_test,windows_FP_test,windows_TN_test,windows_FN_test,num_TP_test, num_FP_test, num_TN_test, num_FN_test,ACC_test,PPV_test,NPV_test,activation_weeks_var1_test,activation_weeks_var2_test] = fun_test(alpha_opt,bg_opt,vec_training_time,windows_test,var1_timeseries,vec_week_test_event)


% Loop on GT data for the selected county. Compute selected windows in
% which there is an increase in the google trends variable
count1 = 0;
count_th_act = 0;
count_deriv_act = 0;

for j1=1:length(windows_test)
    aux0 = var1_timeseries(windows_test{j1})>=bg_opt*max(var1_timeseries(vec_training_time)) ; % pick up gt time series on the chosen window
    aux1 = strfind(aux0',[0 1]) ;                                                                   % indicates how many a crossing (sequnce of [0 1] in aux2) happened: aux3 gives the list where the zeros of [0 1] are located.
    
    aux_windows_test = windows_test{j1};
    alpha_vec_var1  = var1_timeseries(aux_windows_test(2:end))./var1_timeseries(aux_windows_test(1:end-1));
    
    if ~isempty(aux1) && length(find(alpha_vec_var1>1))< alpha_opt                                                    % criteria for event detection: signal should have crossed the threshold at least at 1 point in the sliding window
        count1                             = count1 + 1; % count crossing
        aux_select_windows1                = windows_test{j1};
        select_windows1{count1}            = aux_select_windows1; % save window
        activation_weeks_var1_test{count1} = aux_select_windows1(aux1(1)+1); % save week of first crossing. This is the activation week.
        
        % count threshold activations
        count_th_act = count_th_act +1;
        select_windows_th_act{count_th_act}  = aux_select_windows1; % save window
        
        
        
        clear aux_select_windows1
        
    else if isempty(aux1) && length(find(alpha_vec_var1>1))>=alpha_opt
            count1                        = count1 + 1;
            aux_select_windows1           = windows_test{j1};
            select_windows1{count1}       = aux_select_windows1;                                                 % save window
            
            aux_alpha_vec_var1           = alpha_vec_var1>1;
            aux_ind_alpha_vec            = find(aux_alpha_vec_var1==1);
            aux_ind_activation_alpha     = aux_ind_alpha_vec(alpha_opt);
            activation_weeks_var1_test{count1} = aux_select_windows1(aux_ind_activation_alpha+1);
            
          
            % count derivative activations
            count_deriv_act                         = count_deriv_act +1;
            select_windows_deriv_act{count_deriv_act}  = aux_select_windows1; % save window
            
        else if ~isempty(aux1) && length(find(alpha_vec_var1>1))>=alpha_opt
                
                count1                        = count1 + 1;                                                          % count crossing
                aux_select_windows1           = windows_test{j1};                                                % store selection window
                select_windows1{count1}       = aux_select_windows1;                                                 % save window
                activation_weeks_var1_cross   = aux_select_windows1(aux1(1)+1);                                     % save week of . This is the activation week.
                
                aux_alpha_vec_var1           = alpha_vec_var1>1;
                aux_ind_alpha_vec            = find(aux_alpha_vec_var1==1);
                aux_ind_activation_alpha     = aux_ind_alpha_vec(alpha_opt);
                activation_weeks_var1_alpha  = aux_select_windows1(aux_ind_activation_alpha+1);
                
                activation_weeks_var1_test{count1} = min(activation_weeks_var1_cross,activation_weeks_var1_alpha);
                
            % count threshold activations    
            count_th_act = count_th_act +1;
            select_windows_th_act{count_th_act}  = aux_select_windows1; % save window
        
            
            % count derivative activations
            count_deriv_act                            = count_deriv_act +1;
            select_windows_deriv_act{count_deriv_act}  = aux_select_windows1; % save window
                
            end
        end
    end
    clear aux0 aux1 aux_windows_test alpha_vec_var1
    
end



if count_th_act==0
    select_windows_th_act = cell.empty;
end

if count_deriv_act==0
    select_windows_deriv_act = cell.empty;
end



%% Store events on Rt data for test period.
window_gap = length(windows_test{1})-1;
for j0=1:length(vec_week_test_event)
    
    test_week_event = vec_week_test_event(j0);
    
    if test_week_event>1
        
        aux_count2 = sort(1:min(window_gap,test_week_event-1),'descend');
        
        for count2 = 1:min(window_gap,test_week_event-1)
            aux_select_windows2{j0,aux_count2(count2)} = max(test_week_event-count2,1):test_week_event + window_gap - count2;
            aux_activation_weeks_var2_test{j0,count2} = test_week_event;
        end
    end
    clear test_week_event
end

select_windows2            = reshape(aux_select_windows2',[1,j0*count2]);
activation_weeks_var2_test = reshape(aux_activation_weeks_var2_test',[1,j0*count2]);

%% COMPUTE  OUT OF SAMPLE ACC, PPV AND NPV

% Store info for confusion matrix
num_TP  = 0;
num_TN  = 0;
num_FN  = 0;
num_FP  = 0;


if count1==0
    select_windows1 = cell.empty;
    activation_weeks_var1_test{1} = NaN;
end


% Start Loop on test windows
for j3=1:length(windows_test)
    j3;
    if  ~any(cellfun(@(x) isequal(x, windows_test{j3}), select_windows1)) &&  ~any(cellfun(@(x) isequal(x, windows_test{j3}), select_windows2)) % if there is no var1 AND no var2 activation at the given window, count as true negative
        num_TN             = num_TN+1;
        windows_TN_test{num_TN} = windows_test{j3};
    else if ~any(cellfun(@(x) isequal(x, windows_test{j3}), select_windows1)) &&  any(cellfun(@(x) isequal(x, windows_test{j3}), select_windows2)) % if there is no Var1 activation but var2 activates, counts as false negative
            num_FN = num_FN  + 1;
            windows_FN_test{num_FN} = windows_test{j3};
        else if any(cellfun(@(x) isequal(x, windows_test{j3}), select_windows1)) &&  ~any(cellfun(@(x) isequal(x, windows_test{j3}), select_windows2)) % if there is no Var2 activation but var1 activates, counts as false positive
                num_FP = num_FP  + 1;
                windows_FP_test{num_FP} = windows_test{j3};
                
            else if any(cellfun(@(x) isequal(x, windows_test{j3}), select_windows1)) &&  any(cellfun(@(x) isequal(x, windows_test{j3}), select_windows2)) % if both Var1 and Var2 activates, we ensure that var1 activates BEFORE var2 to count as True positive. Otherwise, we count as false negative. We continue as follows (se lines below)
                    
                    % Start Loop on selected windows for var1
                    for ii1 = 1:length(select_windows1)
                        if isequal(windows_test{j3},select_windows1{ii1}) % check which if the window is exactly a selected window of index ii1
                            dummy_ind1(ii1) = 1;
                        else
                            dummy_ind1(ii1) = 0;
                        end
                    end
                    
                    act_week1= activation_weeks_var1_test{dummy_ind1==1}; % find the activation week.
                    
                    % Start Loop on selected windows for var2
                    for ii2 = 1:length(select_windows2)
                        if isequal(windows_test{j3},select_windows2{ii2}) % check which if the window is exactly a selected window of index ii2
                            dummy_ind2(ii2) = 1;
                        else
                            dummy_ind2(ii2) = 0;
                        end
                    end
                    
                    act_week2  = activation_weeks_var2_test{dummy_ind2==1};% find the activation week.
                    
                    % Check for TP or FN.
                    if act_week1 <  act_week2
                        num_TP = num_TP+1;
                        windows_TP_test{num_TP} = windows_test{j3};
                    else
                        num_FN = num_FN+1;
                        windows_FN_test{num_FN} = windows_test{j3};
                    end
                    
                end
            end
        end
    end
    clear dummy_ind1 dummy_ind2 aux_var1 aux_var2  act_week1  act_week2
end


if num_TP==0
    windows_TP_test = cell.empty;
end


if num_FP==0
    windows_FP_test = cell.empty;
end


if num_TN==0
    windows_TN_test = cell.empty;
end


if num_FN==0
    windows_FN_test = cell.empty;
end


ACC_test = (num_TP + num_TN)/(num_TP + num_TN + num_FP + num_FN); % accuracy
PPV_test = num_TP/(num_TP+num_FP);                                % sensitivity
NPV_test = num_TN/(num_TN + num_FN);                              % specificity

num_TP_test = num_TP;
num_FP_test = num_FP;
num_TN_test = num_TN;
num_FN_test = num_FN;


end


